clc;
clear;

vector_h = [1/2];
a = 0; b = 1000;
% 参数
s = 10; r = 0.03; mu_T = 0.02; mu_b = 0.24; mu_v = 2.4;
k1 = 2.4e-5; k2 = 3e-3; Tmax = 1500;
% 分数阶参数 0.7, 0.8, 0.9, 0.955
alpha = 0.955;

% 二次参数计算
ga = r^alpha/Tmax;
p = r^alpha - mu_T^alpha;
k3 = mu_T^alpha + k2^alpha;
T0 = (p+sqrt(4*ga*s^alpha + p^2))/(2*ga);
Nc = k3*(k1^alpha*T0 + mu_v^alpha)/(k1^alpha*k2^alpha*T0);
ceil(Nc)
N = 6*Nc;

% 初值
tu0 = T0; tl0 = 0; ta0 = 0; v0 = 1e-3;

tic;
error_tu = zeros(length(vector_h),1);
error_tl = zeros(length(vector_h),1);
error_ta = zeros(length(vector_h),1);
error_v = zeros(length(vector_h),1);
for i = 1:length(vector_h)
    h = vector_h(i);
    n = (b-a)/h;
    tn = (a:h:b)';
    % 计算
    [tu, tl, ta, v] = process_Caputo(tu0, tl0, ta0, v0, h, n, s, mu_T, mu_b, mu_v, k1, k2, r, Tmax, N, alpha);
    % [tu, tl, ta, v] = process_ABC(tu0, tl0, ta0, v0, h, n, s, mu_T, mu_b, mu_v, k1, k2, r, Tmax, N, alpha);
    % [tu, tl, ta, v] = process_CF(tu0, tl0, ta0, v0, h, n, s, mu_T, mu_b, mu_v, k1, k2, r, Tmax, N, alpha);
    % 绘图
    process_plot(tn, tu, tl, ta, v);
    % 误差
    if i > 1
        error_tu(i) = process_error(tu_old, tu);
        error_tl(i) = process_error(tl_old, tl);
        error_ta(i) = process_error(ta_old, ta);
        error_v(i) = process_error(v_old, v);
    end
    tu_old = tu;
    tl_old = tl;
    ta_old = ta;
    v_old = v;
    fprintf("已完成\t%d/%d\n", i, length(vector_h));
end
% 误差阶
process_print(error_tu, vector_h);
process_print(error_tl, vector_h);
process_print(error_ta, vector_h);
process_print(error_v, vector_h);

elapsedTOCTime = toc;
disp(["TOC time(s)", num2str(elapsedTOCTime)]);




%% 调用的函数
function [tu, tl, ta, v] = process_Caputo(tu0, tl0, ta0, v0, h, n, s, mu_T, mu_b, mu_v, k1, k2, r, Tmax, N, alpha)
tu = zeros(n+1,1); tl = zeros(n+1,1); ta = zeros(n+1,1); v = zeros(n+1,1);
tu(1) = tu0;
tl(1) = tl0;
ta(1) = ta0;
v(1) = v0;
for k = 1:n
    sum1 = 0;
    sum2 = 0;
    sum3 = 0;
    sum4 = 0;
    for i = 1:k
        sum1 = sum1 + omega(k,i,alpha)*(s^alpha - mu_T^alpha*tu(i) + r^alpha*tu(i)*(1 - (tu(i) + tl(i) + ta(i))/Tmax) - k1^alpha*tu(i)*v(i));
        sum2 = sum2 + omega(k,i,alpha)*(-mu_T^alpha*tl(i) - k2^alpha*tl(i) + k1^alpha*tu(i)*v(i));
        sum3 = sum3 + omega(k,i,alpha)*(-mu_b^alpha*ta(i) + k2^alpha*tl(i));
        sum4 = sum4 + omega(k,i,alpha)*(-mu_v^alpha*v(i) + N*mu_b^alpha*ta(i) - k1^alpha*tu(i)*v(i));
    end
    tu(k+1) = tu(1) + h^alpha/gamma(alpha+1)*sum1;
    tl(k+1) = tl(1) + h^alpha/gamma(alpha+1)*sum2;
    ta(k+1) = ta(1) + h^alpha/gamma(alpha+1)*sum3;
    v(k+1) = v(1) + h^alpha/gamma(alpha+1)*sum4;
end
end

% function [tu, tl, ta, v] = process_CF(tu0, tl0, ta0, v0, h, n, s, mu_T, mu_b, mu_v, k1, k2, r, Tmax, N, alpha)
% tu = zeros(n+1,1); tl = zeros(n+1,1); ta = zeros(n+1,1); v = zeros(n+1,1);
% tu(1) = tu0;
% tl(1) = tl0;
% ta(1) = ta0;
% v(1) = v0;
% for k = 1:n
%     sum1 = 0;
%     sum2 = 0;
%     sum3 = 0;
%     sum4 = 0;
%     for i = 1:k
%         sum1 = sum1 + (s^alpha - mu_T^alpha*tu(i) + r^alpha*tu(i)*(1 - (tu(i) + tl(i) + ta(i))/Tmax) - k1^alpha*tu(i)*v(i));
%         sum2 = sum2 + (-mu_T^alpha*tl(i) - k2^alpha*tl(i) + k1^alpha*tu(i)*v(i));
%         sum3 = sum3 + (-mu_b^alpha*ta(i) + k2^alpha*tl(i));
%         sum4 = sum4 + (-mu_v^alpha*v(i) + N*mu_b^alpha*ta(i) - k1^alpha*tu(i)*v(i));
%     end
%     sum1 = alpha*h*sum1;
%     sum2 = alpha*h*sum2;
%     sum3 = alpha*h*sum3;
%     sum4 = alpha*h*sum4;
%     tu(k+1) = tu(1) + (1-alpha)*() + sum1;
%     tl(k+1) = tl(1) + (1-alpha)*() + sum1;
%     ta(k+1) = ta(1) + (1-alpha)*(-mu_b^alpha*ta(k+1) + k2^alpha*tl(k+1)) + sum1;
%     v(k+1) = v(1) + (1-alpha)*() + sum1;
% end
% end

% function [tu, tl, ta, v] = process_ABC(tu0, tl0, ta0, v0, h, n, s, mu_T, mu_b, mu_v, k1, k2, r, Tmax, N, alpha)
%     tu = zeros(n+1,1); tl = zeros(n+1,1); ta = zeros(n+1,1); v = zeros(n+1,1);
%     tu(1) = tu0;
%     tl(1) = tl0;
%     ta(1) = ta0;
%     v(1) = v0;
%     for k = 1:n
%         sum1 = 0;
%         sum2 = 0;
%         sum3 = 0;
%         sum4 = 0;
%         for i = 1:k
%             sum1 = sum1 + omega(k,i,alpha)*(s^alpha - mu_T^alpha*tu(i) + r^alpha*tu(i)*(1 - (tu(i) + tl(i) + ta(i))/Tmax) - k1^alpha*tu(i)*v(i));
%             sum2 = sum2 + (-mu_T^alpha*tl(i) - k2^alpha*tl(i) + k1^alpha*tu(i)*v(i));
%             sum3 = sum3 + (-mu_b^alpha*ta(i) + k2^alpha*tl(i));
%             sum4 = sum4 + (-mu_v^alpha*v(i) + N*mu_b^alpha*ta(i) - k1^alpha*tu(i)*v(i));
%         end
%         sum1 = (h^alpha)/(beta(alpha)*gamma(alpha))*sum1;
%         sum2 = alpha*h*sum2;
%         sum3 = alpha*h*sum3;
%         sum4 = alpha*h*sum4;
%         tu(k+1) = tu(1) + (1-alpha)/beta(alpha)*() + sum1;
%         tl(k+1) = tl(1) + (1-alpha)*() + sum1;
%         ta(k+1) = ta(1) + (1-alpha)*() + sum1;
%         v(k+1) = v(1) + (1-alpha)*() + sum1;
%     end
% end

function r = omega(k, i, alpha)
r = (k+1-i)^alpha - (k-i)^alpha;
end

function process_plot(tn, tu, tl, ta, v)
figure;
plot(tn, tu, 'r-', 'LineWidth', 1.5); hold on;
plot(tn, tl, 'b--', 'LineWidth', 1.5);
plot(tn, ta, 'k-.', 'LineWidth', 1.5);
plot(tn, v, 'g:', 'LineWidth', 1.5);
xlabel('时间 t');
ylabel('数值解');
legend('tu', 'tl', 'ta', 'v');
title('分数阶模型的数值解');
grid on;
end

function r = process_error(old, new)
r = norm(old - new(1:2:end), inf);
end

function process_print(vector_error, vector_h)
error_var_name = inputname(1);
convergence_rate = zeros(size(vector_h,2), 1);
for i = 2:size(vector_h,2)
    convergence_rate(i) = log(vector_error(i-1)/vector_error(i))/log(vector_h(i-1)/vector_h(i));
end
T = table(compose("1/%1.0f", (1./vector_h)'), compose("%1.4e", vector_error), compose("%1.5f", convergence_rate),...
    'VariableNames', {'h', error_var_name, 'Rate'});
disp(T);
end